1 /* Copyright 2002-2016 CS Systèmes d'Information
2 * Licensed to CS Systèmes d'Information (CS) under one or more
3 * contributor license agreements. See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * CS licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License. You may obtain a copy of the License at
8 *
9 * http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17 package org.orekit.forces.gravity.potential;
18
19 import java.io.IOException;
20 import java.io.InputStream;
21 import java.text.ParseException;
22 import java.util.ArrayList;
23 import java.util.Arrays;
24 import java.util.List;
25 import java.util.Locale;
26
27 import org.hipparchus.util.FastMath;
28 import org.hipparchus.util.Precision;
29 import org.orekit.data.DataLoader;
30 import org.orekit.errors.OrekitException;
31 import org.orekit.errors.OrekitMessages;
32
33 /**This abstract class represents a Gravitational Potential Coefficients file reader.
34 *
35 * <p> As it exits many different coefficients models and containers this
36 * interface represents all the methods that should be implemented by a reader.
37 * The proper way to use this interface is to call the {@link GravityFieldFactory}
38 * which will determine which reader to use with the selected potential
39 * coefficients file.<p>
40 *
41 * @see GravityFieldFactory
42 * @author Fabien Maussion
43 */
44 public abstract class PotentialCoefficientsReader implements DataLoader {
45
46 /** Maximal degree to parse. */
47 private int maxParseDegree;
48
49 /** Maximal order to parse. */
50 private int maxParseOrder;
51
52 /** Regular expression for supported files names. */
53 private final String supportedNames;
54
55 /** Allow missing coefficients in the input data. */
56 private final boolean missingCoefficientsAllowed;
57
58 /** Indicator for complete read. */
59 private boolean readComplete;
60
61 /** Central body reference radius. */
62 private double ae;
63
64 /** Central body attraction coefficient. */
65 private double mu;
66
67 /** Raw tesseral-sectorial coefficients matrix. */
68 private double[][] rawC;
69
70 /** Raw tesseral-sectorial coefficients matrix. */
71 private double[][] rawS;
72
73 /** Indicator for normalized raw coefficients. */
74 private boolean normalized;
75
76 /** Tide system. */
77 private TideSystem tideSystem;
78
79 /** Simple constructor.
80 * <p>Build an uninitialized reader.</p>
81 * @param supportedNames regular expression for supported files names
82 * @param missingCoefficientsAllowed allow missing coefficients in the input data
83 */
84 protected PotentialCoefficientsReader(final String supportedNames,
85 final boolean missingCoefficientsAllowed) {
86 this.supportedNames = supportedNames;
87 this.missingCoefficientsAllowed = missingCoefficientsAllowed;
88 this.maxParseDegree = Integer.MAX_VALUE;
89 this.maxParseOrder = Integer.MAX_VALUE;
90 this.readComplete = false;
91 this.ae = Double.NaN;
92 this.mu = Double.NaN;
93 this.rawC = null;
94 this.rawS = null;
95 this.normalized = false;
96 this.tideSystem = TideSystem.UNKNOWN;
97 }
98
99 /** Get the regular expression for supported files names.
100 * @return regular expression for supported files names
101 */
102 public String getSupportedNames() {
103 return supportedNames;
104 }
105
106 /** Check if missing coefficients are allowed in the input data.
107 * @return true if missing coefficients are allowed in the input data
108 */
109 public boolean missingCoefficientsAllowed() {
110 return missingCoefficientsAllowed;
111 }
112
113 /** Set the degree limit for the next file parsing.
114 * @param maxParseDegree maximal degree to parse (may be safely
115 * set to {@link Integer#MAX_VALUE} to parse all available coefficients)
116 * @since 6.0
117 */
118 public void setMaxParseDegree(final int maxParseDegree) {
119 this.maxParseDegree = maxParseDegree;
120 }
121
122 /** Get the degree limit for the next file parsing.
123 * @return degree limit for the next file parsing
124 * @since 6.0
125 */
126 public int getMaxParseDegree() {
127 return maxParseDegree;
128 }
129
130 /** Set the order limit for the next file parsing.
131 * @param maxParseOrder maximal order to parse (may be safely
132 * set to {@link Integer#MAX_VALUE} to parse all available coefficients)
133 * @since 6.0
134 */
135 public void setMaxParseOrder(final int maxParseOrder) {
136 this.maxParseOrder = maxParseOrder;
137 }
138
139 /** Get the order limit for the next file parsing.
140 * @return order limit for the next file parsing
141 * @since 6.0
142 */
143 public int getMaxParseOrder() {
144 return maxParseOrder;
145 }
146
147 /** {@inheritDoc} */
148 public boolean stillAcceptsData() {
149 return !(readComplete &&
150 getMaxAvailableDegree() >= getMaxParseDegree() &&
151 getMaxAvailableOrder() >= getMaxParseOrder());
152 }
153
154 /** Set the indicator for completed read.
155 * @param readComplete if true, a gravity field has been completely read
156 */
157 protected void setReadComplete(final boolean readComplete) {
158 this.readComplete = readComplete;
159 }
160
161 /** Set the central body reference radius.
162 * @param ae central body reference radius
163 */
164 protected void setAe(final double ae) {
165 this.ae = ae;
166 }
167
168 /** Get the central body reference radius.
169 * @return central body reference radius
170 */
171 protected double getAe() {
172 return ae;
173 }
174
175 /** Set the central body attraction coefficient.
176 * @param mu central body attraction coefficient
177 */
178 protected void setMu(final double mu) {
179 this.mu = mu;
180 }
181
182 /** Get the central body attraction coefficient.
183 * @return central body attraction coefficient
184 */
185 protected double getMu() {
186 return mu;
187 }
188
189 /** Set the {@link TideSystem} used in the gravity field.
190 * @param tideSystem tide system used in the gravity field
191 */
192 protected void setTideSystem(final TideSystem tideSystem) {
193 this.tideSystem = tideSystem;
194 }
195
196 /** Get the {@link TideSystem} used in the gravity field.
197 * @return tide system used in the gravity field
198 */
199 protected TideSystem getTideSystem() {
200 return tideSystem;
201 }
202
203 /** Set the tesseral-sectorial coefficients matrix.
204 * @param rawNormalized if true, raw coefficients are normalized
205 * @param c raw tesseral-sectorial coefficients matrix
206 * (a reference to the array will be stored)
207 * @param s raw tesseral-sectorial coefficients matrix
208 * (a reference to the array will be stored)
209 * @param name name of the file (or zip entry)
210 * @exception OrekitException if a coefficient is missing
211 */
212 protected void setRawCoefficients(final boolean rawNormalized,
213 final double[][] c, final double[][] s,
214 final String name)
215 throws OrekitException {
216
217 // normalization indicator
218 normalized = rawNormalized;
219
220 // set known constant values, if they were not defined in the file.
221 // See Hofmann-Wellenhof and Moritz, "Physical Geodesy",
222 // section 2.6 Harmonics of Lower Degree.
223 // All S_i,0 are irrelevant because they are multiplied by zero.
224 // C0,0 is 1, the central part, since all coefficients are normalized by GM.
225 setIfUnset(c, 0, 0, 1);
226 setIfUnset(s, 0, 0, 0);
227 // C1,0, C1,1, and S1,1 are the x,y,z coordinates of the center of mass,
228 // which are 0 since all coefficients are given in an Earth centered frame
229 setIfUnset(c, 1, 0, 0);
230 setIfUnset(s, 1, 0, 0);
231 setIfUnset(c, 1, 1, 0);
232 setIfUnset(s, 1, 1, 0);
233
234 // cosine part
235 for (int i = 0; i < c.length; ++i) {
236 for (int j = 0; j < c[i].length; ++j) {
237 if (Double.isNaN(c[i][j])) {
238 throw new OrekitException(OrekitMessages.MISSING_GRAVITY_FIELD_COEFFICIENT_IN_FILE,
239 'C', i, j, name);
240 }
241 }
242 }
243 rawC = c;
244
245 // sine part
246 for (int i = 0; i < s.length; ++i) {
247 for (int j = 0; j < s[i].length; ++j) {
248 if (Double.isNaN(s[i][j])) {
249 throw new OrekitException(OrekitMessages.MISSING_GRAVITY_FIELD_COEFFICIENT_IN_FILE,
250 'S', i, j, name);
251 }
252 }
253 }
254 rawS = s;
255
256 }
257
258 /**
259 * Set a coefficient if it has not been set already.
260 * <p>
261 * If {@code array[i][j]} is 0 or NaN this method sets it to {@code value} and returns
262 * {@code true}. Otherwise the original value of {@code array[i][j]} is preserved and
263 * this method return {@code false}.
264 * <p>
265 * If {@code array[i][j]} does not exist then this method returns {@code false}.
266 *
267 * @param array the coefficient array.
268 * @param i degree, the first index to {@code array}.
269 * @param j order, the second index to {@code array}.
270 * @param value the new value to set.
271 * @return {@code true} if the coefficient was set to {@code value}, {@code false} if
272 * the coefficient was not set to {@code value}. A {@code false} return indicates the
273 * coefficient has previously been set to a non-NaN, non-zero value.
274 */
275 private boolean setIfUnset(final double[][] array,
276 final int i,
277 final int j,
278 final double value) {
279 if (array.length > i && array[i].length > j &&
280 (Double.isNaN(array[i][j]) || Precision.equals(array[i][j], 0.0, 0))) {
281 // the coefficient was not already initialized
282 array[i][j] = value;
283 return true;
284 } else {
285 return false;
286 }
287 }
288
289 /** Get the maximal degree available in the last file parsed.
290 * @return maximal degree available in the last file parsed
291 * @since 6.0
292 */
293 public int getMaxAvailableDegree() {
294 return rawC.length - 1;
295 }
296
297 /** Get the maximal order available in the last file parsed.
298 * @return maximal order available in the last file parsed
299 * @since 6.0
300 */
301 public int getMaxAvailableOrder() {
302 return rawC[rawC.length - 1].length - 1;
303 }
304
305 /** {@inheritDoc} */
306 public abstract void loadData(InputStream input, String name)
307 throws IOException, ParseException, OrekitException;
308
309 /** Get a provider for read spherical harmonics coefficients.
310 * @param wantNormalized if true, the provider will provide normalized coefficients,
311 * otherwise it will provide un-normalized coefficients
312 * @param degree maximal degree
313 * @param order maximal order
314 * @return a new provider
315 * @exception OrekitException if the requested maximal degree or order exceeds the
316 * available degree or order or if no gravity field has read yet
317 * @see #getConstantProvider(boolean, int, int)
318 * @since 6.0
319 */
320 public abstract RawSphericalHarmonicsProvider getProvider(boolean wantNormalized, int degree, int order)
321 throws OrekitException;
322
323 /** Get a time-independent provider for read spherical harmonics coefficients.
324 * @param wantNormalized if true, the raw provider must provide normalized coefficients,
325 * otherwise it will provide un-normalized coefficients
326 * @param degree maximal degree
327 * @param order maximal order
328 * @return a new provider, with no time-dependent parts
329 * @exception OrekitException if the requested maximal degree or order exceeds the
330 * available degree or order or if no gravity field has read yet
331 * @see #getProvider(boolean, int, int)
332 * @since 6.0
333 */
334 protected ConstantSphericalHarmonics getConstantProvider(final boolean wantNormalized,
335 final int degree, final int order)
336 throws OrekitException {
337
338 if (!readComplete) {
339 throw new OrekitException(OrekitMessages.NO_GRAVITY_FIELD_DATA_LOADED);
340 }
341
342 if (degree >= rawC.length) {
343 throw new OrekitException(OrekitMessages.TOO_LARGE_DEGREE_FOR_GRAVITY_FIELD,
344 degree, rawC.length - 1);
345 }
346
347 if (order >= rawC[rawC.length - 1].length) {
348 throw new OrekitException(OrekitMessages.TOO_LARGE_ORDER_FOR_GRAVITY_FIELD,
349 order, rawC[rawC.length - 1].length);
350 }
351
352 // fix normalization
353 final double[][] truncatedC = buildTriangularArray(degree, order, 0.0);
354 final double[][] truncatedS = buildTriangularArray(degree, order, 0.0);
355 rescale(1.0, normalized, rawC, rawS, wantNormalized, truncatedC, truncatedS);
356
357 return new ConstantSphericalHarmonics(ae, mu, tideSystem, truncatedC, truncatedS);
358
359 }
360
361 /** Build a coefficients triangular array.
362 * @param degree array degree
363 * @param order array order
364 * @param value initial value to put in array elements
365 * @return built array
366 */
367 protected static double[][] buildTriangularArray(final int degree, final int order, final double value) {
368 final int rows = degree + 1;
369 final double[][] array = new double[rows][];
370 for (int k = 0; k < array.length; ++k) {
371 array[k] = buildRow(k, order, value);
372 }
373 return array;
374 }
375
376 /**
377 * Parse a double from a string. Accept the Fortran convention of using a 'D' or
378 * 'd' instead of an 'E' or 'e'.
379 *
380 * @param string to be parsed.
381 * @return the double value of {@code string}.
382 */
383 protected static double parseDouble(final String string) {
384 return Double.parseDouble(string.toUpperCase(Locale.ENGLISH).replace('D', 'E'));
385 }
386
387 /** Build a coefficients row.
388 * @param degree row degree
389 * @param order row order
390 * @param value initial value to put in array elements
391 * @return built row
392 */
393 protected static double[] buildRow(final int degree, final int order, final double value) {
394 final double[] row = new double[FastMath.min(order, degree) + 1];
395 Arrays.fill(row, value);
396 return row;
397 }
398
399 /** Extend a list of lists of coefficients if needed.
400 * @param list list of lists of coefficients
401 * @param degree degree required to be present
402 * @param order order required to be present
403 * @param value initial value to put in list elements
404 */
405 protected void extendListOfLists(final List<List<Double>> list, final int degree, final int order,
406 final double value) {
407 for (int i = list.size(); i <= degree; ++i) {
408 list.add(new ArrayList<Double>());
409 }
410 final List<Double> listN = list.get(degree);
411 final Double v = Double.valueOf(value);
412 for (int j = listN.size(); j <= order; ++j) {
413 listN.add(v);
414 }
415 }
416
417 /** Convert a list of list into an array.
418 * @param list list of lists of coefficients
419 * @return a new array
420 */
421 protected double[][] toArray(final List<List<Double>> list) {
422 final double[][] array = new double[list.size()][];
423 for (int i = 0; i < array.length; ++i) {
424 array[i] = new double[list.get(i).size()];
425 for (int j = 0; j < array[i].length; ++j) {
426 array[i][j] = list.get(i).get(j);
427 }
428 }
429 return array;
430 }
431
432 /** Parse a coefficient.
433 * @param field text field to parse
434 * @param list list where to put the coefficient
435 * @param i first index in the list
436 * @param j second index in the list
437 * @param cName name of the coefficient
438 * @param name name of the file
439 * @exception OrekitException if the coefficient is already set
440 */
441 protected void parseCoefficient(final String field, final List<List<Double>> list,
442 final int i, final int j,
443 final String cName, final String name)
444 throws OrekitException {
445 final double value = parseDouble(field);
446 final double oldValue = list.get(i).get(j);
447 if (Double.isNaN(oldValue) || Precision.equals(oldValue, 0.0, 0)) {
448 // the coefficient was not already initialized
449 list.get(i).set(j, value);
450 } else {
451 throw new OrekitException(OrekitMessages.DUPLICATED_GRAVITY_FIELD_COEFFICIENT_IN_FILE,
452 name, i, j, name);
453 }
454 }
455
456 /** Parse a coefficient.
457 * @param field text field to parse
458 * @param array array where to put the coefficient
459 * @param i first index in the list
460 * @param j second index in the list
461 * @param cName name of the coefficient
462 * @param name name of the file
463 * @exception OrekitException if the coefficient is already set
464 */
465 protected void parseCoefficient(final String field, final double[][] array,
466 final int i, final int j,
467 final String cName, final String name)
468 throws OrekitException {
469 final double value = parseDouble(field);
470 final double oldValue = array[i][j];
471 if (Double.isNaN(oldValue) || Precision.equals(oldValue, 0.0, 0)) {
472 // the coefficient was not already initialized
473 array[i][j] = value;
474 } else {
475 throw new OrekitException(OrekitMessages.DUPLICATED_GRAVITY_FIELD_COEFFICIENT_IN_FILE,
476 name, i, j, name);
477 }
478 }
479
480 /** Rescale coefficients arrays.
481 * @param scale general scaling factor to apply to all elements
482 * @param normalizedOrigin if true, the origin coefficients are normalized
483 * @param originC cosine part of the origina coefficients
484 * @param originS sine part of the origin coefficients
485 * @param wantNormalized if true, the rescaled coefficients must be normalized
486 * @param rescaledC cosine part of the rescaled coefficients to fill in (may be the originC array)
487 * @param rescaledS sine part of the rescaled coefficients to fill in (may be the originS array)
488 * @exception OrekitException if normalization/unnormalization fails because of an underflow
489 * due to too high degree/order
490 */
491 protected static void rescale(final double scale,
492 final boolean normalizedOrigin, final double[][] originC,
493 final double[][] originS, final boolean wantNormalized,
494 final double[][] rescaledC, final double[][] rescaledS)
495 throws OrekitException {
496
497 if (wantNormalized == normalizedOrigin) {
498 // apply only the general scaling factor
499 for (int i = 0; i < rescaledC.length; ++i) {
500 final double[] rCi = rescaledC[i];
501 final double[] rSi = rescaledS[i];
502 final double[] oCi = originC[i];
503 final double[] oSi = originS[i];
504 for (int j = 0; j < rCi.length; ++j) {
505 rCi[j] = oCi[j] * scale;
506 rSi[j] = oSi[j] * scale;
507 }
508 }
509 } else {
510
511 // we have to re-scale the coefficients
512 // (we use rescaledC.length - 1 for the order instead of rescaledC[rescaledC.length - 1].length - 1
513 // because typically trend or pulsation arrays are irregular, some test cases have
514 // order 2 elements at degree 2, but only order 1 elements for higher degrees for example)
515 final double[][] factors = GravityFieldFactory.getUnnormalizationFactors(rescaledC.length - 1,
516 rescaledC.length - 1);
517
518 if (wantNormalized) {
519 // normalize the coefficients
520 for (int i = 0; i < rescaledC.length; ++i) {
521 final double[] rCi = rescaledC[i];
522 final double[] rSi = rescaledS[i];
523 final double[] oCi = originC[i];
524 final double[] oSi = originS[i];
525 final double[] fi = factors[i];
526 for (int j = 0; j < rCi.length; ++j) {
527 final double factor = scale / fi[j];
528 rCi[j] = oCi[j] * factor;
529 rSi[j] = oSi[j] * factor;
530 }
531 }
532 } else {
533 // un-normalize the coefficients
534 for (int i = 0; i < rescaledC.length; ++i) {
535 final double[] rCi = rescaledC[i];
536 final double[] rSi = rescaledS[i];
537 final double[] oCi = originC[i];
538 final double[] oSi = originS[i];
539 final double[] fi = factors[i];
540 for (int j = 0; j < rCi.length; ++j) {
541 final double factor = scale * fi[j];
542 rCi[j] = oCi[j] * factor;
543 rSi[j] = oSi[j] * factor;
544 }
545 }
546 }
547
548 }
549 }
550
551 }